H = H_Int,
v = v,
zeta = zeta,
pup = 0.5)^2),
lower = p0 * (1 - pup),
upper = p0 * (1 + pup))
P_Int <- I$par
S_Int <- share_2(p = P_Int, deltar = DATA$deltar, v = v, zeta = zeta)
# --------------------------------------------------------------------------
# SECCIÓN 6: SIMULACIÓN DE LOS EFECTOS COORDINADOS (CARTEL - ESCENARIO TIGO)
# --------------------------------------------------------------------------
# Propósito: Simular un escenario de colusión/cartelización incorporando a
# COMCEL junto con las firmas integradas (TIGO y MOVISTAR).
# Se reasignan las identidades de las firmas para reflejar el cartel
dx <- DATA %>%
mutate(firm_id = ifelse(EMPRESA %in% c("COMCEL", "TIGO", "MOVISTAR"), "CARTEL", EMPRESA))
# Creamos nuevamente la matriz de propiedad H para el escenario del cartel
H_coor <- matrix(0, nrow = nrow(DATA), ncol = nrow(DATA))
for(m_idx in seq_along(unique(t))) {
market_rows <- which(mkti == m_idx)
# Extraemos las firmas de este mercado
firms_in_market <- dx$firm_id[market_rows]
matching_firms <- outer(firms_in_market, firms_in_market, "==") * 1
H_coor[market_rows, market_rows] <- matching_firms
}
# Ejecutar la simulación de efectos coordinados
p0 <- p
pup <- 0.5
C <- nlminb(start = P_Int,
objective = function(x) sum(mktsimul(x, p0 = p0,
deltar = DATA$deltar,
cm = DATA$cm,
H = H_coor,
v = v,
zeta = zeta,
pup = 0.5)^2),
lower = p0 * (1 - pup),
upper = p0 * (1 + pup))
P_coor <- C$par
S_coor <- share_2(p = P_coor, deltar = DATA$deltar, v = v, zeta = zeta)
# --------------------------------------------------------------------------
# SECCIÓN 7.1: ALTERNATIVA REGULATORIA 1
# --------------------------------------------------------------------------
# Propósito: Simular la situación bajo la Alternativa 1
PA1 <- 2.5
x <- expand.grid(
p = PA1,
FECHA = ts,
EMPRESA = unique(d4$EMPRESA)
) %>%
left_join(
distinct(d4, EMPRESA, AS),
by = "EMPRESA",
relationship = "many-to-many"
)
x$L_T <- predict(R1, newdata = x)
x$T_A <- exp(x$L_T)
d <- filter(d4, FECHA %in% ts)
d <- select(d, FECHA, EMPRESA, AS, T)
x <- merge(x, d, by = c("FECHA", "EMPRESA", "AS"), all.x=TRUE)
x <- x %>% filter(!is.na(T))
x$d <- x$T_A - x$T
y <- x %>%
group_by(FECHA, EMPRESA) %>%
summarise(
D = sum(d, na.rm = TRUE),
.groups = 'drop'
) %>%
rename(t = FECHA)
d <- DATA1
d <- merge(d, y, by = c("t", "EMPRESA"), all.x = TRUE)
d <- d %>% mutate(D = coalesce(D, 0), DD = D + TA2P)
d <- d %>% mutate(A2P = DD/q)
d <- d %>% mutate(A2P = ifelse(A2P < 0, 0, A2P))
d <- d %>%
mutate(
EMPRESA = recode(EMPRESA, "MOVISTAR" = "TIGO")
) %>%
mutate(s = S_Int) %>%
group_by(t, EMPRESA, posp) %>%
mutate(
A2P = if_else(EMPRESA == "TIGO",
weighted.mean(A2P, s, na.rm = TRUE),
A2P)
) %>%
ungroup()
XB <- XPF[i,]
XB[,3] <- d$A2P
deltarA1 <- as.vector(XB %*% b)
xx <- data.frame(EMPRESA = d$EMPRESA, posp = d$posp, deltar = deltarA1, s = S_Int)
xx <- xx %>%
group_by(EMPRESA, posp) %>%
mutate(
deltar = if_else(EMPRESA == "TIGO",
weighted.mean(deltar, s, na.rm = TRUE),
deltar)
) %>%
ungroup()
deltarA1 <- as.vector(xx$deltar)
# Usamos la matriz de propiedad H como una matriz de identidad por bloques H_Int
p0 <- p
pup <- 0.5
IA1 <- nlminb(start = p,
objective = function(x) sum(mktsimul(x, p0 = p0,
deltar = deltarA1,
cm = DATA$cm,
H = H_Int,
v = v,
zeta = zeta,
pup = 0.5)^2),
lower = p0 * (1 - pup),
upper = p0 * (1 + pup))
P_IntA1 <- IA1$par
S_IntA1 <- share_2(p = P_IntA1, deltar = deltarA1, v = v, zeta = zeta)
# --------------------------------------------------------------------------
# SECCIÓN 7.2: ALTERNATIVA REGULATORIA 1 CON CARTEL
# --------------------------------------------------------------------------
# Propósito: Simular un escenario de colusión/cartelización incorporando a
# COMCEL junto con las firmas integradas (TIGO y MOVISTAR)
# Ejecutar la simulación de efectos coordinados utilizando H_coor como matriz
# de propiedad en el cartel
p0 <- p
pup <- 0.5
CA1 <- nlminb(start = P_IntA1,
objective = function(x) sum(mktsimul(x, p0 = p0,
deltar = deltarA1,
cm = DATA$cm,
H = H_coor,
v = v,
zeta = zeta,
pup = 0.5)^2),
lower = p0 * (1 - pup),
upper = p0 * (1 + pup))
P_coorA1 <- CA1$par
S_coorA1 <- share_2(p = P_coorA1, deltar = deltarA1, v = v, zeta = zeta)
# --------------------------------------------------------------------------
# SECCIÓN 8.1: ALTERNATIVA REGULATORIA 2
# --------------------------------------------------------------------------
# Propósito: Simular la situación bajo la Alternativa 2
PA2 <- 8.098
x <- expand.grid(
p = PA2,
FECHA = ts,
EMPRESA = unique(d4$EMPRESA)
) %>%
left_join(
distinct(d4, EMPRESA, AS),
by = "EMPRESA",
relationship = "many-to-many"
)
x$L_T <- predict(R1, newdata = x)
x$T_A <- exp(x$L_T)
d <- filter(d4, FECHA %in% ts)
d <- select(d, FECHA, EMPRESA, AS, T)
x <- merge(x, d, by = c("FECHA", "EMPRESA", "AS"), all.x=TRUE)
x <- x %>% filter(!is.na(T))
x$d <- x$T_A - x$T
y <- x %>%
group_by(FECHA, EMPRESA) %>%
summarise(
D = sum(d, na.rm = TRUE),
.groups = 'drop'
) %>%
rename(t = FECHA)
d <- DATA1
d <- merge(d, y, by = c("t", "EMPRESA"), all.x = TRUE)
d <- d %>% mutate(D = coalesce(D, 0), DD = D + TA2P)
d <- d %>% mutate(A2P = DD/q)
d <- d %>% mutate(A2P = ifelse(A2P < 0, 0, A2P))
d <- d %>%
mutate(
EMPRESA = recode(EMPRESA, "MOVISTAR" = "TIGO")
) %>%
group_by(t, EMPRESA, posp) %>%
mutate(
A2P = if_else(EMPRESA == "TIGO",
weighted.mean(A2P, s, na.rm = TRUE),
A2P)
) %>%
ungroup() %>%
XB <- XPF[i,]
EMPRESA <- DATA$EMPRESA
DATA <- DATA %>% arrange(t, EMPRESA, serv, posp)
X <- DATA %>%
transmute(
posp = posp, # Dummy de pospago
A2P = A2P
) %>% as.matrix()
# --- Generación de Efectos Fijos por empresa, año y fecha ---
dx <- DATA %>%
transmute(
p,
XF = paste0('_', EMPRESA), # Efectos fijos por firma
XT = paste0('_', ANNO)    # Efectos fijos por año
)
# Generación de matrices binarias eliminando la primera columna para evitar multicolinealidad perfecta ([, -1])
XF <- model.matrix(~ 0 + XF, data = dx)[,-1]
XT <- model.matrix(~ 0 + XT, data = dx)[,-1]
# --------------------------------------------------------------------------
# SECCIÓN 0: FUNCIONES DE PARTICIPACIONES PREDICHAS
# --------------------------------------------------------------------------
# Propósito: Calcular las cuotas de mercado predichas.
## Función auxiliar para calcular las participaciones de mercado (shares)
share <- function(delta, mui){
num   <- exp(delta + mui) # Numerador de la probabilidad Logit
# Calcular el denominador (suma inclusiva) por mercado
den   <- t(sapply(unique(t), function(m) colSums(num[which(t == m), ])))[mkti, ]
si    <- num / (1 + den)  # Probabilidades individuales
s     <- rowMeans(si)     # Cuotas de mercado agregadas predichas
return(list(si = si, s = s))
}
share_2 <- function(p, deltar, v, zeta){
## Parámetro de precio
alpha <- zeta[2]
## Efectos de precio
mui   <- - zeta[1]*v*p
delta <- deltar - alpha*p
## Recálculo de cuotas
s <- share(delta = delta, mui = mui)
si  <- s$si
s   <- s$s
return(s)
}
# --------------------------------------------------------------------------
# SECCIÓN 1: FUNCIONES DE PUNTO FIJO (CONTRACCIÓN DE BERRY)
# --------------------------------------------------------------------------
# Propósito: Calcular el vector de utilidades medias (delta) que iguala las
# cuotas de mercado observadas con las predichas por el modelo BLP.
## Función de punto fijo para delta
# delta: vector de utilidades medias actuales
# mui: parte no lineal de la utilidad (interacción consumidor-producto)
fixedn <- function(delta, mui) {
# Participaciones predichas
s <- share(delta = delta, mui = mui)
s <- s$s
# Actualización de la contracción de Berry
delta <- as.numeric(delta + log(S) - log(s))
return(delta)
}
## Función de punto fijo para zeta (parámetros no lineales)
# zeta: vector de parámetros no lineales a estimar
# Utiliza el algoritmo 'squarem' para acelerar la convergencia del punto fijo.
fixed <- function(zeta) {
# Construcción de 'mui' (utilidad idiosincrásica) basada en características individuales
mui   <- zeta[1]*v*p
# Intento de encontrar delta usando squarem
delta <- try(squarem(delta0, fixedn, mui = mui, control = list(tol = 1e-10)), silent = T)
expl  <- class(delta)
if(expl == 'try-error') {
si <- delta <- NULL # Fallo en la convergencia
} else {
delta <- delta$par  # delta convergido
si <- share(delta = delta, mui = mui)
si <- si$si
}
return(list(delta = delta, si = si, expl = expl))
}
# --------------------------------------------------------------------------
# SECCIÓN 2: RECUPERACIÓN DE PRIMITIVAS DEL MODELO
# --------------------------------------------------------------------------
# Propósito: Dado un vector de parámetros zeta, recuperar los parámetros lineales
# (beta, gamma), costos marginales y errores estructurales (xi, omega).
struct <- function(zeta) {
## Punto fijo para delta
Fout  <- fixed(zeta)
expl  <- Fout$expl
if(expl == "try-error") {
out <- list(expl = expl)
} else {
## Delta y cuotas
delta  <- Fout$delta
si     <- Fout$si        				  # probabilidades individuales
s      <- rowMeans(si)  				  # cuotas de mercado
## Parámetros lineales de demanda (IV GMM)
alpha  <- zeta[2]
deltar <- delta + alpha*p
b <- XZPZ %*% deltar                                  # Estimador IV para beta
## Márgenes y Costos:
# Elasticidad precio derivada
aim         <- -(alpha + zeta[1]*v)
dsdp        <- -(aim * si) %*% t(si) / NS             # Derivada cruzada ds_r/dp_j
diag(dsdp)  <- rowMeans(aim * si * (1 - si))          # Derivada propia ds_j/dp_j
# Inversión de la matriz de derivadas para obtener márgenes (Nash-Bertrand)
TH          <- H * dsdp                               # H es la matriz de propiedad (ownership matrix)
mkup        <- -as.numeric(solve(TH, tol = 1e-40) %*% s)
cm          <- p - mkup                               # Costos marginales recuperados
## Parámetros lineales de oferta (IV GMM)
g <- WZPZ %*% cm         				  # Estimador IV para gamma (costos)
## Residuales estructurales
xi <- deltar - XPF %*% b 				  # Error de demanda
wo <- cm - WPF %*% g     				  # Error de oferta (costos)
## Salida: parámetros y primitivas
out <- list(b = b, g = g, delta = delta, cm = cm, si = si, xi = xi, wo = wo, expl = expl)
}
return(out)
}
# --------------------------------------------------------------------------
# SECCIÓN 3: ESTIMACIÓN GMM (MÉTODO GENERALIZADO DE MOMENTOS)
# --------------------------------------------------------------------------
# Propósito: Función objetivo a minimizar para encontrar los parámetros óptimos.
## Función de momentos GMM
gmm <- function(zeta) {
## Restricciones:
print(zeta) 			   # Imprimir iteración actual para monitoreo
## Obtener primitivas
strc <- struct(zeta)
expl <- strc$expl
if(expl == 'try-error') {
out <- 1e+10 		   # Penalización grande si falla el punto fijo
} else {
## Residuales
xi <- strc$xi
cm <- strc$cm
wo <- strc$wo
# Penalizaciones adicionales para costos marginales ilógicos (negativos o > precio)
wo <- wo + ((cm > p)*(cm - p)^2 + (cm < 0)*(100*cm)^2)
## Valor de la función objetivo
om  <- c(xi, wo) 		   # Vector apilado de errores
out <- t(om) %*% ZPZ2 %*% om   # Forma cuadrática GMM (m' W m)
out <- as.numeric(out)/1e+5    # Escalado
}
print(out)
return(out)
}
## Restricciones sobre los parámetros
cons <- function(zeta) {
zeta     <- abs(zeta) # Forzar positividad (ej. para desviaciones estándar)
return(zeta)
}
# --------------------------------------------------------------------------
# SECCIÓN 4: FUNCIONES DE ERROR PARA ERRORES ESTÁNDAR
# --------------------------------------------------------------------------
# Propósito: Funciones auxiliares para calcular la varianza de los estimadores
# (probablemente usadas para calcular el Jacobiano numérico).
## Funciones de error respecto a theta_1 (parámetros lineales demanda)
omega1 <- function(beta) {
alpha   <- zeta[2]
delta_r <- delta + alpha*p
xi <- delta_r - XPF %*% beta
wo <- wo
return(c(xi, wo))
}
## Funciones de error respecto a theta_2 (parámetros no lineales)
omega2 <- function(zeta) {
out  <- struct(zeta)
xi   <- out$xi
wo   <- out$wo
return(c(xi, wo))
}
# --------------------------------------------------------------------------
# SECCIÓN 5: SIMULACIÓN DE MERCADO (CONTRAFACTUALES)
# --------------------------------------------------------------------------
# Propósito: Calcular nuevos precios de equilibrio ante cambios en el mercado
# (ej. fusión) iterando sobre las condiciones de primer orden (FOC).
mktsimul <- function(p, p0, deltar, cm, H, v, zeta, pup) {
## Restricciones sobre los valores de p
# p <- pmin(p, (1 + pup) * p0)
# p <- pmax(p, (1 - pup) * p0)
# p <- pmax(p, 10)
## Parámetros de precio
alpha <- zeta[2]
aim   <- -(alpha + zeta[1]*v)
## Efectos de precio
mui   <- - zeta[1]*v*p
delta <- deltar - alpha*p
## Recálculo de cuotas
s <- share(delta = delta, mui = mui)
si  <- s$si
s   <- s$s
## Recálculo de márgenes y precios (FOC de oferta):
dsdp        <- -(aim * si) %*% t(si) / NS
diag(dsdp)  <- rowMeans(aim * si * (1 - si))
TH          <- H * dsdp
mkup        <- -as.numeric(solve(TH, tol = 1e-40) %*% s)
eq          <- p - cm - mkup
## Salida
return(eq)
}
p_simul <- function(p, pup){
p <- pmin(p, (1 + pup) * p0)
p <- pmax(p, (1 - pup) * p0)
p <- pmax(p, 10)
return(p)
}
# --------------------------------------------------------------------------
# SECCIÓN 6: CÁLCULO DE MEJOR RESPUESTA A PARTIR DE LAS FOC
# --------------------------------------------------------------------------
best_response_2 <- function(p_optim) {
# Aumentar el contador cada vez que optim llama a la función
iter_count <<- iter_count + 1
p <- numeric(nrow(a))
p[idx_change] <- p_optim
p[idx_fixed]  <- p_fixed
mui    <- - zeta[1]*v*p
deltar <- a$deltar
delta  <- deltar - zeta[2]*p
## Recálculo de cuotas
s <- share(delta = delta, mui = mui)
si  <- s$si
s   <- s$s
s_optim <- s[idx_change]
M <- a$M[idx_change]
cm <- a$CM[idx_change]
pi_optim <- (p_optim - cm) %*% (M*s_optim)
cat(sprintf("Iter: %d | Beneficio: %.4f | Precio Promedio: %.2f\n",
iter_count, pi_optim, mean(p_optim)))
if(!is.finite(pi_optim)) return(-1e10)
return(as.numeric(pi_optim))
}
# --------------------------------------------------------------------------
# SECCIÓN 7: CÁLCULO DE BIENESTAR
# --------------------------------------------------------------------------
welfare <- function(p, deltar, v, zeta, M){
aim <- zeta[2] + zeta[1]*v
aim <- aim[1:length(ts), ] # Ajuste de dimensiones si es necesario
mui <- -zeta[1]*v*p
delta  <- deltar - zeta[2]*p
# Cálculo del valor esperado de la utilidad máxima (Log-sum / Inclusive Value)
num <- exp(delta + mui)
x   <- t(sapply(unique(t), function(m) colSums(num[which(t == m), ])))
# Bienestar del consumidor (Consumer Surplus)
W  <- sum(M*rowMeans(log(x)/aim, na.rm = T))/1000
return(W)
}
library(readxl)     # Para lectura de archivos Excel
library(dplyr)      # Para manipulación y agregación de datos
library(lubridate)  # Para manipulación de fechas
library(writexl)    # Para exportar resultados a Excel
load("C:/Users/maria.marino/OneDrive - COMISION DE REGULACIÓN DE COMUNICACIONES/CRC/Fraude/Bases/Script a/c. v.RData")
load("C:/Users/maria.marino/OneDrive - COMISION DE REGULACIÓN DE COMUNICACIONES/CRC/Fraude/Bases/Script a/c. Estimout_01.RData")
load("C:/Users/maria.marino/COMISION DE REGULACIÓN DE COMUNICACIONES/[2025] [DR] Medidas para mitigar el fraude cibernético - Documentos/2. Documento Amarillo/17. Documentos Fase 2/3. Anexos documento soporte/Anexo Evaluación de alternativas_Publico/a. BD.RData")
load("C:/Users/maria.marino/COMISION DE REGULACIÓN DE COMUNICACIONES/[2025] [DR] Medidas para mitigar el fraude cibernético - Documentos/2. Documento Amarillo/17. Documentos Fase 2/3. Anexos documento soporte/ANEXOE~2/AREGRE~1.RDA")
load("C:/Users/maria.marino/COMISION DE REGULACIÓN DE COMUNICACIONES/[2025] [DR] Medidas para mitigar el fraude cibernético - Documentos/2. Documento Amarillo/17. Documentos Fase 2/3. Anexos documento soporte/Anexo Evaluación de alternativas_Publico/b. MDATA4.RData")
load("C:/Users/maria.marino/OneDrive - COMISION DE REGULACIÓN DE COMUNICACIONES/CRC/Fraude/Bases/Script a/c. v.RData")
load("C:/Users/maria.marino/OneDrive - COMISION DE REGULACIÓN DE COMUNICACIONES/CRC/Fraude/Bases/Script a/c. Estimout_01.RData")
load("C:/Users/maria.marino/COMISION DE REGULACIÓN DE COMUNICACIONES/[2025] [DR] Medidas para mitigar el fraude cibernético - Documentos/2. Documento Amarillo/17. Documentos Fase 2/3. Anexos documento soporte/Anexo Evaluación de alternativas_Publico/a. BD.RData")
load("C:/Users/maria.marino/COMISION DE REGULACIÓN DE COMUNICACIONES/[2025] [DR] Medidas para mitigar el fraude cibernético - Documentos/2. Documento Amarillo/17. Documentos Fase 2/3. Anexos documento soporte/Anexo Evaluación de alternativas_Publico/b. MDATA4.RData")
load("C:/Users/maria.marino/COMISION DE REGULACIÓN DE COMUNICACIONES/[2025] [DR] Medidas para mitigar el fraude cibernético - Documentos/2. Documento Amarillo/17. Documentos Fase 2/3. Anexos documento soporte/ANEXOE~2/AREGRE~1.RDA")
EMPRESA <- DATA$EMPRESA
DATA <- DATA %>% arrange(t, EMPRESA, serv, posp)
X <- DATA %>%
transmute(
posp = posp, # Dummy de pospago
A2P = A2P
) %>% as.matrix()
# --- Generación de Efectos Fijos por empresa, año y fecha ---
dx <- DATA %>%
transmute(
p,
XF = paste0('_', EMPRESA), # Efectos fijos por firma
XT = paste0('_', ANNO)    # Efectos fijos por año
)
# Generación de matrices binarias eliminando la primera columna para evitar multicolinealidad perfecta ([, -1])
XF <- model.matrix(~ 0 + XF, data = dx)[,-1]
XT <- model.matrix(~ 0 + XT, data = dx)[,-1]
XPF <- cbind(constant = 1, X, XF, XT)
DATA1 <- DATA
# 1. Filtrado y transformación inicial
# Aprovechamos para filtrar DATA y sincronizar v, delta y cm en un solo paso
i <- which(DATA$t %in% ts)
EMPRESA <- DATA$EMPRESA
DATA <- DATA %>% arrange(t, EMPRESA, serv, posp)
X <- DATA %>%
transmute(
posp = posp, # Dummy de pospago
A2P = A2P
) %>% as.matrix()
# --- Generación de Efectos Fijos por empresa, año y fecha ---
dx <- DATA %>%
transmute(
p,
XF = paste0('_', EMPRESA), # Efectos fijos por firma
XT = paste0('_', ANNO)    # Efectos fijos por año
)
# Generación de matrices binarias eliminando la primera columna para evitar multicolinealidad perfecta ([, -1])
XF <- model.matrix(~ 0 + XF, data = dx)[,-1]
XT <- model.matrix(~ 0 + XT, data = dx)[,-1]
XPF <- cbind(constant = 1, X, XF, XT)
DATA1 <- DATA
# --------------------------------------------------------------------------
# SECCIÓN 3: PREPARACIÓN DE DATOS
# --------------------------------------------------------------------------
# Propósito: Filtrar la base para el periodo de simulación, sincronizar los
# vectores (v, delta, cm) y extraer los parámetros base.
# 1. Filtrado y transformación inicial
# Aprovechamos para filtrar DATA y sincronizar v, delta y cm en un solo paso
i <- which(DATA$t %in% ts)
library(readxl)     # Para lectura de archivos Excel
library(dplyr)      # Para manipulación y agregación de datos
library(lubridate)  # Para manipulación de fechas
library(writexl)    # Para exportar resultados a Excel
DATA <- DATA %>% arrange(t, EMPRESA, serv, posp)
X <- DATA %>%
transmute(
posp = posp, # Dummy de pospago
A2P = A2P
) %>% as.matrix()
# --- Generación de Efectos Fijos por empresa, año y fecha ---
dx <- DATA %>%
transmute(
p,
XF = paste0('_', EMPRESA), # Efectos fijos por firma
XT = paste0('_', ANNO)    # Efectos fijos por año
)
# Generación de matrices binarias eliminando la primera columna para evitar multicolinealidad perfecta ([, -1])
XF <- model.matrix(~ 0 + XF, data = dx)[,-1]
XT <- model.matrix(~ 0 + XT, data = dx)[,-1]
XPF <- cbind(constant = 1, X, XF, XT)
DATA1 <- DATA
# 1. Filtrado y transformación inicial
# Aprovechamos para filtrar DATA y sincronizar v, delta y cm en un solo paso
i <- which(DATA$t %in% ts)
